
//Adds the first robust first stage f-stat as described in Andrews, Stock, and Sun (2019) p.737; stores result in "fStat"
program addRobustFStat
	version 16.1
	args varToAdd
	estadd scalar fStat = (_b[`varToAdd']^2)/e(V)["`varToAdd'", "`varToAdd'"]
	
	
end

//Adds the first treatment coefficient and standard error; stores result in "treatb" and "treatse"
program addTreat
	version 16.1
	args varToAdd
	qui estadd local treatse = "(0"+string(round(_se[`varToAdd'], .001))+")"
	
	
	local tempb = string(round(_b[`varToAdd'], .001))
	
	di "hello"
	
	local colonIndex = ustrrpos("`varToAdd'", ":")
	local strLen = strlen("`varToAdd'")
	local eqno = udsubstr("`varToAdd'",1,`colonIndex'-1)
	local varName = udsubstr("`varToAdd'",`colonIndex'+1,`strLen'+`colonIndex'-1)
	di "`colonIndex'"
	di "`eqno'"
	di "`varName'"
	
	test [`eqno']`varName'
	di "hello2"
	
	
	//Determine Stars for alternate null
	local stars = ""
			
	if(`r(p)'<0.01){
		local stars = "***"
	}	
	else if(`r(p)'<0.05){
		local stars = "**"
	}
	else if(`r(p)'<0.1){
		local stars = "*"
	}

	estadd local treatb = "0`tempb'"+"`stars'"
	
	
end

//Program to test some specified null; 
program testNull, rclass
	version 16.1
	args pi delta piVar deltaVar covar null
	
	tempname arStat pVal
	
	scalar `arStat' = ((`delta'-`pi'*`null')^2)/(`deltaVar'+(`null')^2*`piVar'-2*`null'*`covar')
	scalar `pVal' = 1-chi2(1, `arStat')
	
	
	return scalar pi = `pi'
	return scalar delta = `delta'
	return scalar piVar = `piVar'
	return scalar deltaVar = `deltaVar'
	return scalar covar = `covar'
	return scalar arStat = `arStat'
	return scalar pVal = `pVal'
	return scalar null = `null'
	return scalar beta = `delta'/`pi'

end

//Program to jointly estimate first stage and reduced form; argument structure is first-stage model, reduced form model, cluster variable, treatmentVar, and the model title
program jointEstimate, rclass
	version 16.1 //Establish Version
	//Parse arguments
	gettoken fsVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken clusterVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatmentVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken modelTitle 0:0, parse(",")

	di "Working on:`modelTitle'"
	
	
	
	//Establish Temporary Variables
	tempname fsK rfK N nClust fsCorrection rfCorrection pi delta piVarUncorrected deltaVarUncorrected covarUncorrected piVar deltaVar covar covarHigh covarLow
	
	regress `fsVars'
	scalar `fsK' = e(df_m)+1
	return scalar fsK = `fsK'
	estimates store fs 
		
	regress `rfVars'
	scalar `rfK' = e(df_m)+1
	return scalar rfK = `rfK'
	estimates store rf 
	
	qui suest fs rf, vce(cluster `clusterVar')
	scalar `nClust' = e(N_clust)
	scalar `N' = e(N)
	scalar `pi' = e(b)["y1", "fs_mean:`treatmentVar'"]
	scalar `delta' = e(b)["y1", "rf_mean:`treatmentVar'"]
	scalar `piVarUncorrected' = e(V)["fs_mean:`treatmentVar'", "fs_mean:`treatmentVar'"]
	scalar `deltaVarUncorrected' = e(V)["rf_mean:`treatmentVar'", "rf_mean:`treatmentVar'"]
	scalar `covarUncorrected' = e(V)["rf_mean:`treatmentVar'", "fs_mean:`treatmentVar'"]
	
	
	return scalar nClust = `nClust'
	return scalar N = `N'
	scalar `fsCorrection' = (`N'-1)/(`N'-`fsK'+(`nClust'-1))
	scalar `rfCorrection' = (`N'-1)/(`N'-`rfK'+(`nClust'-1))
	scalar `piVar' = `piVarUncorrected'*`fsCorrection'
	scalar `deltaVar' = `deltaVarUncorrected'*`fsCorrection'
	scalar `covarHigh' = max(`fsCorrection'*`covarUncorrected',`rfCorrection'*`covarUncorrected')
	scalar `covarLow' = min(`fsCorrection'*`covarUncorrected',`rfCorrection'*`covarUncorrected')
	scalar `covar' = `covarLow'
	
	return scalar fsCorrection = `fsCorrection'
	return scalar rfCorrection = `rfCorrection'
	return scalar pi = `pi'
	return scalar delta = `delta'
	return scalar piVarUncorrected = `piVarUncorrected'
	return scalar deltaVarUncorrected = `deltaVarUncorrected'
	return scalar covarUncorrected = `covarUncorrected'
	return scalar piVar = `piVar'
	return scalar deltaVar = `deltaVar'
	return scalar covarHigh = `covarHigh'
	return scalar covarLow = `covarLow'
	return scalar covar = `covar'
	
	
end

program regFWL, eclass
	version 16.1 //Establish Version
	
	gettoken out_var 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken est_var 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken partial_var 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken absorb_vars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken cluster_var 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken res_var 0:0, parse(",")
	gettoken com 0:0, parse(",")
	
	
	if("`out_var'"=="_") local out_var = ""
	if("`est_var'"=="_") local est_var = ""
	if("`partial_var'"=="_") local partial_var = ""
	if("`absorb_vars'"=="_") local absorb_vars = ""
	if("`cluster_var'"=="_") local cluster_var = ""
	if("`res_var'"=="_") local res_var = ""
	
	local numClusterVar = 0
	
	foreach tempVar in `cluster_var'{
		local numClusterVar = `numClusterVar'+1
		local clusterVar`numClusterVar' = "`tempVar'"
	}
	di "absorb: `absorb_vars'"
	qui reghdfe `out_var' `partial_var', absorb(`absorb_vars') res(`out_var'_res)
	local new_dep = "`out_var'_res"
	local new_ind = ""
	
	local num_est = 0
	if("`est_var'"!=""){
		foreach temp_var in `est_var'{
			qui reghdfe `temp_var' `partial_var', absorb(`absorb_vars') res(`temp_var'_res)
			local new_ind = "`new_ind' `temp_var'_res"
			local num_est = `num_est'+1
		}
	}
	
	//Run just to get correct values for k; this is the full regression, so other statistics that are unrelated to covariance can be used
	if("`res_var'"==""){
		qui reghdfe `out_var' `est_var' `partial_var', absorb(`absorb_vars') cluster(`cluster_var')
	}
	else{
		qui reghdfe `out_var' `est_var' `partial_var', absorb(`absorb_vars') cluster(`cluster_var') res(`res_var')
	}
	local k = e(df_a)+e(df_m)
	local N = e(N)
	local r2 = e(r2)
	qui reghdfe `new_dep' `new_ind', noabsorb nocon cluster(`cluster_var')
	
	qui regress `new_dep' `new_ind', nocon
	ereturn list
	
	estimates store mod1
	
	qui suest mod1, vce(cluster `clusterVar1')
	local nClust1 = e(N_clust)
	mat var1 = e(V)
	
	qui suest mod1, vce(cluster `clusterVar2')
	local nClust2 = e(N_clust)
	mat var2 = e(V)
	
	qui suest mod1, vce(robust)
	
	mat varR = e(V)
	mat beta = e(b)
	
	
	mat mat_var = var1+var2-varR
	local nClust = min(`nClust1', `nClust2')
	local correction = (`nClust')/(`nClust'-1)*(`N'-1)/(`N'-`k')
	mat mat_var = `correction'*mat_var
	matlist mat_var
	ereturn post beta mat_var
	ereturn scalar r2 = `r2'
	ereturn scalar N = `N'
	ereturn scalar N_clust = `nClust'
	drop `new_dep' `new_ind'

	//Complete
end



program jointEstimateFWL, rclass
	version 16.1 //Establish Version
	//Parse arguments
	gettoken fsOut 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfOut 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken fsControls 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfControls 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatControl 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken absorbVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken clusterVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatmentVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken modelTitle 0:0, parse(",")

	di "Working on:`modelTitle'"
	
	local numClusterVar = 0
	
	foreach tempVar in `clusterVar'{
		local numClusterVar = `numClusterVar'+1
		local clusterVar`numClusterVar' = "`tempVar'"
	}
	
	//Establish Temporary Variables
	tempname fsK rfK N nClust nClust1 nClust2 fsCorrection rfCorrection pi delta piVarUncorrected deltaVarUncorrected covarUncorrected piVar deltaVar covar covarHigh covarLow piVarUncorrected1 deltaVarUncorrected1 covarUncorrected1 piVarUncorrected2 deltaVarUncorrected2 covarUncorrected2 piVarUncorrectedR deltaVarUncorrectedR covarUncorrectedR 
	
	qui reghdfe `fsOut' `fsControls', absorb(`absorbVars') res(`fsOut'_res) cluster(`clusterVar')
	scalar `fsK' = e(df_a)+e(df_m)+1 //Adding 1 because second regression only has one regressor by construction, so this accounts for all regressors
	
	qui reghdfe `treatControl' `fsControls', absorb(`absorbVars') res(fs_`treatControl'_res)
	
	regress `fsOut'_res fs_`treatControl'_res, nocon
	
	return scalar fsK = `fsK'
	estimates store fs 
		
	qui reghdfe `rfOut' `rfControls', absorb(`absorbVars') res(`rfOut'_res) cluster(`clusterVar')
	scalar `rfK' = e(df_a)+e(df_m)+1 //Adding 1 because second regression only has one regressor by construction, so this accounts for all regressors
	
	qui reghdfe `treatControl' `rfControls', absorb(`absorbVars') res(rf_`treatControl'_res)
	
	qui regress `rfOut'_res rf_`treatControl'_res, nocon
	
	
	return scalar rfK = `rfK'
	estimates store rf 
	
	if `numClusterVar'==1{
		qui suest fs rf, vce(cluster `clusterVar')
	
		scalar `nClust' = e(N_clust)
		scalar `N' = e(N)
		scalar `pi' = e(b)["y1", "fs_mean:fs_`treatControl'_res"]
		scalar `delta' = e(b)["y1", "rf_mean:rf_`treatControl'_res"]
		scalar `piVarUncorrected' = e(V)["fs_mean:fs_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		scalar `deltaVarUncorrected' = e(V)["rf_mean:rf_`treatControl'_res", "rf_mean:rf_`treatControl'_res"]
		scalar `covarUncorrected' = e(V)["rf_mean:rf_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
	}
	else{
		di "Current approach can handle two-way clustering. Higher dimensional clustering will produce incorrect results."
		qui suest fs rf, vce(cluster `clusterVar1')
		scalar `nClust1' = e(N_clust)
		scalar `N' = e(N)
		scalar `pi' = e(b)["y1", "fs_mean:fs_`treatControl'_res"]
		scalar `delta' = e(b)["y1", "rf_mean:rf_`treatControl'_res"]
		scalar `piVarUncorrected1' = e(V)["fs_mean:fs_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		scalar `deltaVarUncorrected1' = e(V)["rf_mean:rf_`treatControl'_res", "rf_mean:rf_`treatControl'_res"]
		scalar `covarUncorrected1' = e(V)["rf_mean:rf_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		
		return scalar piVarUncorrected1 = `piVarUncorrected1'
		return scalar deltaVarUncorrected1 = `deltaVarUncorrected1'
		return scalar covarUncorrected1 = `covarUncorrected1'
		
		qui suest fs rf, vce(cluster `clusterVar2')
		
		scalar `nClust2' = e(N_clust)
		scalar `piVarUncorrected2' = e(V)["fs_mean:fs_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		scalar `deltaVarUncorrected2' = e(V)["rf_mean:rf_`treatControl'_res", "rf_mean:rf_`treatControl'_res"]
		scalar `covarUncorrected2' = e(V)["rf_mean:rf_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		
		return scalar piVarUncorrected2 = `piVarUncorrected2'
		return scalar deltaVarUncorrected2 = `deltaVarUncorrected2'
		return scalar covarUncorrected2 = `covarUncorrected2'
		
		qui suest fs rf, vce(robust)
		
		scalar `piVarUncorrectedR' = e(V)["fs_mean:fs_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		scalar `deltaVarUncorrectedR' = e(V)["rf_mean:rf_`treatControl'_res", "rf_mean:rf_`treatControl'_res"]
		scalar `covarUncorrectedR' = e(V)["rf_mean:rf_`treatControl'_res", "fs_mean:fs_`treatControl'_res"]
		
		return scalar piVarUncorrectedR = `piVarUncorrectedR'
		return scalar deltaVarUncorrectedR = `deltaVarUncorrectedR'
		return scalar covarUncorrectedR = `covarUncorrectedR'
		
		scalar `nClust' = min(`nClust1', `nClust2')
		return scalar nClust1 = `nClust1'
		return scalar nClust2 = `nClust2'
		scalar `piVarUncorrected' = `piVarUncorrected1'+`piVarUncorrected2'-`piVarUncorrectedR'
		scalar `deltaVarUncorrected' = `deltaVarUncorrected1'+`deltaVarUncorrected2'-`deltaVarUncorrectedR'
		scalar `covarUncorrected' = `covarUncorrected1'+`covarUncorrected2'-`covarUncorrectedR'
	}
	
	
	
	return scalar nClust = `nClust'
	return scalar N = `N'
	scalar `fsCorrection' = (`nClust')/(`nClust'-1)*(`N'-1)/(`N'-`fsK')
	scalar `rfCorrection' = (`nClust')/(`nClust'-1)*(`N'-1)/(`N'-`rfK')
	scalar `piVar' = `piVarUncorrected'*`fsCorrection'
	scalar `deltaVar' = `deltaVarUncorrected'*`fsCorrection'
	scalar `covarHigh' = max(`fsCorrection'*`covarUncorrected',`rfCorrection'*`covarUncorrected')
	scalar `covarLow' = min(`fsCorrection'*`covarUncorrected',`rfCorrection'*`covarUncorrected')
	scalar `covar' = `covarLow'
	
	return scalar fsCorrection = `fsCorrection'
	return scalar rfCorrection = `rfCorrection'
	return scalar pi = `pi'
	return scalar delta = `delta'
	return scalar piVarUncorrected = `piVarUncorrected'
	return scalar deltaVarUncorrected = `deltaVarUncorrected'
	return scalar covarUncorrected = `covarUncorrected'
	return scalar piVar = `piVar'
	return scalar deltaVar = `deltaVar'
	return scalar covarHigh = `covarHigh'
	return scalar covarLow = `covarLow'
	return scalar covar = `covar'
	
	drop `fsOut'_res fs_`treatControl'_res `rfOut'_res rf_`treatControl'_res
	
end


//Program to output jointly estimated resuls along with a range of null hypotheses; argument structure is first-stage model, reduced form model, cluster variable, treatmentVar, modelTitle, nullLow, nullDelta, nullHigh, primary null
program runModel
	version 16.1
	gettoken fsVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken clusterVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatmentVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken modelTitle 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullLow 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullDelta 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullHigh 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken primaryNull 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull2 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull3 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull4 0:0, parse(",")
	
	tempname null tempPi tempDelta tempPiVar tempDeltaVar tempCovar fileName
	
	
	if("`c(os)'"=="Windows"){
		local tableSaveDir = "Tables\"
		local plotSaveDir = "Plots\"
		local pythonSaveDir = "PythonScripts\"
	}
	else{
		local tableSaveDir = "Tables/"
		local plotSaveDir = "Plots/"
		local pythonSaveDir = "PythonScripts/"
	}
	
	local fileName = "`tableSaveDir'iv-"+"`modelTitle'"+".xlsx"
	
	jointEstimate `fsVars', `rfVars', `clusterVar', `treatmentVar', `modelTitle'
	return list
	
	scalar `tempPi' = `r(pi)'
	scalar `tempDelta' = `r(delta)'
	scalar `tempPiVar' = `r(piVar)'
	scalar `tempDeltaVar' = `r(deltaVar)'
	scalar `tempCovar' = `r(covar)'
	
	python script "`pythonSaveDir'calculateARStat.py", args(`tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `nullLow' `nullDelta' `nullHigh' `modelTitle')
	
	putexcel set "`fileName'", sheet("JointResults") modify
	
	local primaryNullString = string(`primaryNull', "%2.1f")
	local altNullString = string(`altNull', "%2.1f")
	
	qui putexcel B1 = `r(pi)'
	qui putexcel B2 = `r(delta)'
	qui putexcel B3 = `r(piVar)'
	qui putexcel B4 = `r(deltaVar)'
	qui putexcel B5 = `r(covar)'
	qui putexcel B6 = `r(covarHigh)'
	qui putexcel B7 = `r(covarLow)'
	qui putexcel B8 = `r(piVarUncorrected)'
	qui putexcel B9 = `r(deltaVarUncorrected)'
	qui putexcel B10 = `r(covarUncorrected)'
	qui putexcel B11 = `r(fsCorrection)'
	qui putexcel B12 = `r(rfCorrection)'
	qui putexcel B13 = `r(N)'
	qui putexcel B14 = `r(nClust)'
	qui putexcel B15 = `r(fsK)'
	qui putexcel B16 = `r(rfK)'
	
	qui putexcel A1 = "pi"
	qui putexcel A2 = "delta"
	qui putexcel A3 = "pi variance"
	qui putexcel A4 = "delta variance"
	qui putexcel A5 = "pi-delta covariance"
	qui putexcel A6 = "pi-delta covariance high"
	qui putexcel A7 = "pi-delta covariance low"
	qui putexcel A8 = "pi variance uncorrected"
	qui putexcel A9 = "delta variance uncorrected"
	qui putexcel A10 = "pi-delta covariance uncorrected"
	qui putexcel A11 = "first stage finite-sample correction"
	qui putexcel A12 = "reduced form finite-sample correction"
	qui putexcel A13 = "N"
	qui putexcel A14 = "N Clusters"
	qui putexcel A15 = "first stage k"
	qui putexcel A16 = "reduced form k"
	qui putexcel A17 = "AR(Beta=`primaryNullString')"
	qui putexcel A18 = "P-Value(Beta=`primaryNullString')"
	qui putexcel A19 = "beta"
	qui putexcel A20 = "AR(Beta=`altNullString')"
	qui putexcel A21 = "P-Value(Beta=`altNullString')"
	qui putexcel A22 = "AR(Beta=`altNull2')"
	qui putexcel A23 = "P-Value(Beta=`altNull2')"
	qui putexcel A24 = "AR(Beta=`altNull3')"
	qui putexcel A25 = "P-Value(Beta=`altNull3')"
	qui putexcel A26 = "AR(Beta=`altNull4')"
	qui putexcel A27 = "P-Value(Beta=`altNull4')"
	
	
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `primaryNull'
	
	qui putexcel B17 = `r(arStat)'
	qui putexcel B18 = `r(pVal)'
	qui putexcel B19 = `r(beta)'
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull'
	qui putexcel B20 = `r(arStat)'
	qui putexcel B21 = `r(pVal)'
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull2'
	qui putexcel B22 = `r(arStat)'
	qui putexcel B23 = `r(pVal)'
	
	if("`altNull2'"!="`altNull3'"){
		qui putexcel A24 = "AR(Beta=`altNull3')"
		qui putexcel A25 = "P-Value(Beta=`altNull3')"
		
		testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull3'
		qui putexcel B24 = `r(arStat)'
		qui putexcel B25 = `r(pVal)'
	}
	if("`altNull3'"!="`altNull4'"){
		qui putexcel A26 = "AR(Beta=`altNull4')"
		qui putexcel A27 = "P-Value(Beta=`altNull4')"
		
		testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull4'
		qui putexcel B26 = `r(arStat)'
		qui putexcel B27 = `r(pVal)'
	}
	

end

//Program to output jointly estimated resuls along with a range of null hypotheses; argument structure is first-stage model, reduced form model, cluster variable, treatmentVar, modelTitle, nullLow, nullDelta, nullHigh, primary null
program runModelFWL
	version 16.1
	gettoken fsOut 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfOut 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken fsControls 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken rfControls 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatControl 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken absorbVars 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken clusterVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken treatmentVar 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken modelTitle 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullLow 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullDelta 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken nullHigh 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken primaryNull 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull2 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull3 0:0, parse(",")
	gettoken com 0:0, parse(",")
	gettoken altNull4 0:0, parse(",")
	
	tempname null tempPi tempDelta tempPiVar tempDeltaVar tempCovar fileName
	
	
	if("`c(os)'"=="Windows"){
		local tableSaveDir = "Tables\"
		local plotSaveDir = "Plots\"
		local pythonSaveDir = "PythonScripts\"
	}
	else{
		local tableSaveDir = "Tables/"
		local plotSaveDir = "Plots/"
		local pythonSaveDir = "PythonScripts/"
	}
	
	local fileName = "`tableSaveDir'iv-"+"`modelTitle'"+".xlsx"
	
	di "`fsOut'"
	di "`rfOut'"
	di "`fsControls'"
	di "`rfControls'"
	di "`treatControl'"
	di "`absorbVars'"
	
	jointEstimateFWL `fsOut', `rfOut', `fsControls', `rfControls', `treatControl', `absorbVars', `clusterVar', `treatmentVar', `modelTitle'
	return list
	
	scalar `tempPi' = `r(pi)'
	scalar `tempDelta' = `r(delta)'
	scalar `tempPiVar' = `r(piVar)'
	scalar `tempDeltaVar' = `r(deltaVar)'
	scalar `tempCovar' = `r(covar)'
	
	python script "`pythonSaveDir'calculateARStat.py", args(`tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `nullLow' `nullDelta' `nullHigh' `modelTitle')
	
	putexcel set "`fileName'", sheet("JointResults") modify
	
	local primaryNullString = string(`primaryNull', "%2.1f")
	local altNullString = string(`altNull', "%2.1f")
	
	qui putexcel B1 = `r(pi)'
	qui putexcel B2 = `r(delta)'
	qui putexcel B3 = `r(piVar)'
	qui putexcel B4 = `r(deltaVar)'
	qui putexcel B5 = `r(covar)'
	qui putexcel B6 = `r(covarHigh)'
	qui putexcel B7 = `r(covarLow)'
	qui putexcel B8 = `r(piVarUncorrected)'
	qui putexcel B9 = `r(deltaVarUncorrected)'
	qui putexcel B10 = `r(covarUncorrected)'
	qui putexcel B11 = `r(fsCorrection)'
	qui putexcel B12 = `r(rfCorrection)'
	qui putexcel B13 = `r(N)'
	qui putexcel B14 = `r(nClust)'
	qui putexcel B15 = `r(fsK)'
	qui putexcel B16 = `r(rfK)'
	
	qui putexcel A1 = "pi"
	qui putexcel A2 = "delta"
	qui putexcel A3 = "pi variance"
	qui putexcel A4 = "delta variance"
	qui putexcel A5 = "pi-delta covariance"
	qui putexcel A6 = "pi-delta covariance high"
	qui putexcel A7 = "pi-delta covariance low"
	qui putexcel A8 = "pi variance uncorrected"
	qui putexcel A9 = "delta variance uncorrected"
	qui putexcel A10 = "pi-delta covariance uncorrected"
	qui putexcel A11 = "first stage finite-sample correction"
	qui putexcel A12 = "reduced form finite-sample correction"
	qui putexcel A13 = "N"
	qui putexcel A14 = "N Clusters"
	qui putexcel A15 = "first stage k"
	qui putexcel A16 = "reduced form k"
	qui putexcel A17 = "AR(Beta=`primaryNullString')"
	qui putexcel A18 = "P-Value(Beta=`primaryNullString')"
	qui putexcel A19 = "beta"
	qui putexcel A20 = "AR(Beta=`altNullString')"
	qui putexcel A21 = "P-Value(Beta=`altNullString')"
	qui putexcel A22 = "AR(Beta=`altNull2')"
	qui putexcel A23 = "P-Value(Beta=`altNull2')"
	qui putexcel A24 = "AR(Beta=`altNull3')"
	qui putexcel A25 = "P-Value(Beta=`altNull3')"
	qui putexcel A26 = "AR(Beta=`altNull4')"
	qui putexcel A27 = "P-Value(Beta=`altNull4')"
	
	
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `primaryNull'
	
	qui putexcel B17 = `r(arStat)'
	qui putexcel B18 = `r(pVal)'
	qui putexcel B19 = `r(beta)'
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull'
	qui putexcel B20 = `r(arStat)'
	qui putexcel B21 = `r(pVal)'
	
	testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull2'
	qui putexcel B22 = `r(arStat)'
	qui putexcel B23 = `r(pVal)'
	
	if("`altNull2'"!="`altNull3'"){
		qui putexcel A24 = "AR(Beta=`altNull3')"
		qui putexcel A25 = "P-Value(Beta=`altNull3')"
		
		testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull3'
		qui putexcel B24 = `r(arStat)'
		qui putexcel B25 = `r(pVal)'
	}
	if("`altNull3'"!="`altNull4'"){
		qui putexcel A26 = "AR(Beta=`altNull4')"
		qui putexcel A27 = "P-Value(Beta=`altNull4')"
		
		testNull `tempPi' `tempDelta' `tempPiVar' `tempDeltaVar' `tempCovar' `altNull4'
		qui putexcel B26 = `r(arStat)'
		qui putexcel B27 = `r(pVal)'
	}
	

end

//Construct table for LaTex output
program constructTable
	version 16.1
	
	gettoken tableLabel 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken outcomeLabel 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken models 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken modelTitles 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken primaryNull 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken altNull 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken altNull2 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken altNull3 0:0, parse(";")
	gettoken com 0:0, parse(";")
	gettoken altNull4 0:0, parse(";")
	

	if("`c(os)'"=="Windows"){
		local tableSaveDir = "Tables\"
		local plotSaveDir = "Plots\"
		local pythonSaveDir = "PythonScripts\"
	}
	else{
		local tableSaveDir = "Tables/"
		local plotSaveDir = "Plots/"
		local pythonSaveDir = "PythonScripts/"
	}
	
	local i = 1
	
	local floatFormat = "%04.3f"
	
	foreach model in `models'{
		di "`model'"
		preserve
		
		clear
		
		qui import excel stat=A arg1=B arg2=C arg3=D arg4=E using "`tableSaveDir'iv.xlsx", sheet("`model'") cellrange(A1:E29)
	
	
		gen y = 1
		gen x = 1
		
		qui reg y x
		qui eststo `model'
		
		estadd local groupFE = "Y"
		estadd local timeFE = "Y"
		estadd local timeTrend = "Y"
		estadd local commonControls = "Y"
		
		
		//Geta Beta point estimate
		qui sum arg1 if stat=="beta"
		qui estadd local betaPoint = string(round(r(mean), .001), "`floatFormat'")
		qui estadd local modelNum = "(`i')"
		
		if(`i'==1){
			qui estadd local dsc = "None"
		} 
		else if(`i'==2){
			qui estadd local dsc = "Documented"
		} 
		else if(`i'==3){
			qui estadd local dsc = "\makecell{Large \nextline Undocumented}"
		} 
		else if(`i'==4){
			qui estadd local dsc = "\makecell{Small \nextline Undocumented}"
		}
		
		local i = `i'+1
		
		local floatFormat = "%04.3f"
		
		//Get P-Value and Stars for third alternate null (if appropriate)
		if("`tableLabel'"== "GasUseiv" || "`tableLabel'"=="GasUsePerCapitaiv"){
			//Get primary null hypothesis test info
			//Get AR-Stat primary null
			
			qui sum arg1 if stat=="P-Value(Beta=11.0)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for primary null
			qui estadd local primaryP = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=11.0)"
			//Determine Stars for primary null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local primaryARStat = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			qui sum arg1 if stat=="P-Value(Beta=9.7)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for primary null
			qui estadd local altP = string(round(r(mean),0.001), "`floatFormat'")
			
			qui sum arg1 if stat=="AR(Beta=9.7)"
			//Determine Stars for primary null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			//Get alternate null hypothesis test info
			//Get AR-Stat alternate null
			qui sum arg1 if stat=="P-Value(Beta=1)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for alternate null
			qui estadd local altP2 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=1)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat2 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			
			//Get P-Value and Stars for second alternate null
			qui sum arg1 if stat=="P-Value(Beta=-4.8)"
			qui local tempP = r(mean)
			qui estadd local altP3 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=-4.8)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat3 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			
			//Get P-Value and Stars for second alternate null
			qui sum arg1 if stat=="P-Value(Beta=-6.1)"
			qui local tempP = r(mean)
			qui estadd local altP4 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=-6.1)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat4 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
		}
		else{
			
			//Get primary null hypothesis test info
			//Get AR-Stat primary null
			
			qui sum arg1 if stat=="P-Value(Beta=1.6)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for primary null
			qui estadd local primaryP = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=1.6)"
			//Determine Stars for primary null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local primaryARStat = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			//Get alternate null hypothesis test info
			//Get AR-Stat alternate null
			qui sum arg1 if stat=="P-Value(Beta=1.0)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for alternate null
			qui estadd local altP = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=1.0)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			
			//Get alternate null hypothesis test info
			//Get AR-Stat alternate null
			qui sum arg1 if stat=="P-Value(Beta=.3)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for alternate null
			qui estadd local altP2 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=.3)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat2 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			
			//Get alternate null hypothesis test info
			//Get AR-Stat alternate null
			qui sum arg1 if stat=="P-Value(Beta=0)"
			qui local tempP = r(mean)
			//Get P-Value and Stars for alternate null
			qui estadd local altP3 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=0)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat3 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
			
			//Get P-Value and Stars for second alternate null
			qui sum arg1 if stat=="P-Value(Beta=-2.5)"
			qui local tempP = r(mean)
			qui estadd local altP4 = string(round(r(mean),0.001), "`floatFormat'")
			qui sum arg1 if stat=="AR(Beta=-2.5)"
			//Determine Stars for alternate null
			qui local stars = ""
			
			if(`tempP'<0.01){
				local stars = "***"
			}	
			else if(`tempP'<0.05){
				local stars = "**"
			}
			else if(`tempP'<0.1){
				local stars = "*"
			}
			qui estadd local altARStat4 = string(round(r(mean),.001), "`floatFormat'")+"`stars'"
		}
		//Get N
		qui sum arg1 if stat=="N"
		qui estadd local obs = string(round(r(mean)), "%9.0fc")
		
		//Get 95% Confidence set
		local floatFormat = "%3.2f"
		qui sum arg3 if stat=="95% Confidence Set"
		if(r(N)==1){
			qui sum arg2 if stat=="95% Confidence Set"
			local lPoint = string(round(r(mean),.01), "`floatFormat'")
			sum arg3 if stat=="95% Confidence Set"
			qui local rPoint = string(round(r(mean),.01), "`floatFormat'")
			
			qui estadd local cs95 "\makecell{(-$\infty$, $`lPoint']$; \nextline $[`rPoint'$, $\infty$)}"
			
		}
		else{
			qui sum arg1 if stat=="95% Confidence Set"
			if(r(N)==1){
				local lPoint = string(round(r(mean),.01), "`floatFormat'")
			}
			else{
				local lPoint = "-$\infty$"
			}
			qui sum arg2 if stat=="95% Confidence Set"
			if(r(N)==1){
				local rPoint = string(round(r(mean),.01), "`floatFormat'")
			}
			else{
				local rPoint = "$\infty$"
			}
			
			if("`lPoint'"=="-$\infty$" & "`rPoint'"=="$\infty$"){
				qui estadd local cs95 = "(`lPoint', `rPoint')"
			}
			else{
				qui estadd local cs95 = "[`lPoint', `rPoint']"
			}
			
			
			
		}
		
		//Get 90% Confidence Set
		qui sum arg3 if stat=="90% Confidence Set"
		if(r(N)==1){
			qui sum arg2 if stat=="90% Confidence Set"
			local lPoint = string(round(r(mean),.01))
			qui sum arg3 if stat=="90% Confidence Set"
			local rPoint = string(round(r(mean),.01))
			
			qui estadd local cs90 "(-$\infty$, `lPoint']; [`rPoint', $\infty$)"
			
		}
		else{
			qui sum arg1 if stat=="90% Confidence Set"
			if(r(N)==1){
				local lPoint = string(round(r(mean),.01))
			}
			else{
				local lPoint = "-$\infty$"
			}
			qui sum arg2 if stat=="90% Confidence Set"
			if(r(N)==1){
				local rPoint = string(round(r(mean),.01))
			}
			else{
				local rPoint = "$\infty$"
			}
			
			qui estadd local cs90 = "[`lPoint', `rPoint']"
			
		}
		
		
		
		
		restore
	}
	
	//Output Table
	
	if("`c(os)'"=="Windows"){
		local tableSaveDir = "Tables\IV\"
		local plotSaveDir = "Plots\"
		local pythonSaveDir = "PythonScripts\"
	}
	else{
		local tableSaveDir = "Tables/IV/"
		local plotSaveDir = "Plots/"
		local pythonSaveDir = "PythonScripts/"
	}
	
	
	if("`tableLabel'"=="VMTiv"){
		local titleOutcome = "Log Registrations on Log VMT"
		local noteString = "\hline \hline \end{tabular} \vspace{0.2cm} \begin{tablenotes} \begin{spacing}{0} \begin{footnotesize} \item \textbf{Notes:} The top panel of this table presents IV estimates and 95\% confidence sets robust to weak instruments of the linear relationship between log vehicle-miles travelled (VMT) and log vehicle registrations.  The bottom panel presents Anderson-Rubin chi-squared statistics along with p-values for different null hypotheses pertaining to the parameter estimated in the top panel.  The unit of observation is state-year.  All regressions include state fixed effects, year fixed effects and state-specific linear time trends. The four specifications pertaining to different levels of ``Data Source Controls'' (DSCs) and control specifications are documented in Table \cref{ssec:Spec}. Confidence sets are based on standard errors clustered by state. * p\textless 0.1, ** p\textless 0.05, *** p\textless 0.01 \end{footnotesize} \end{spacing} \end{tablenotes} \end{table}"
	} 
	else if("`tableLabel'"=="VMTPerCapitaiv"){
		local titleOutcome = "Log Registrations per Capita on Log VMT per Capita"
		local noteString = "\hline \hline \end{tabular} \vspace{0.2cm} \begin{tablenotes} \begin{spacing}{0} \begin{footnotesize} \item \textbf{Notes:} The top panel of this table presents IV estimates and 95\% confidence sets robust to weak instruments of the linear relationship between log vehicle-miles travelled (VMT) per capita and log vehicle registrations per capita.  The bottom panel presents Anderson-Rubin chi-squared statistics along with p-values for different null hypotheses pertaining to the parameter estimated in the top panel.  The unit of observation is state-year. The four specifications pertaining to different levels of ``Data Source Controls'' (DSCs) and control specifications are documented in Table \cref{ssec:Spec}. Confidence sets are based on standard errors clustered by state. * p\textless 0.1, ** p\textless 0.05, *** p\textless 0.01 \end{footnotesize} \end{spacing} \end{tablenotes} \end{table}"
	}
	else if("`tableLabel'"=="GasUseiv"){
		local titleOutcome = "Log Registrations on Log Highway Gas Use"
		local noteString = "\hline \hline \end{tabular} \vspace{0.2cm} \begin{tablenotes} \begin{spacing}{0} \begin{footnotesize} \item \textbf{Notes:} The top panel of this table presents IV estimates and 95\% confidence sets robust to weak instruments of the linear relationship between log highway gas use and log vehicle registrations.  The bottom panel presents Anderson-Rubin chi-squared statistics along with p-values for different null hypotheses pertaining to the parameter estimated in the top panel. The unit of observation is state-year. The four specifications pertaining to different levels of ``Data Source Controls'' (DSCs) and control specifications are documented in Table \cref{ssec:Spec}. Confidence sets are based on standard errors clustered by state. * p\textless 0.1, ** p\textless 0.05, *** p\textless 0.01 \end{footnotesize} \end{spacing} \end{tablenotes} \end{table}"
	}
	else if("`tableLabel'"=="GasUsePerCapitaiv"){
		local titleOutcome = "Log Registrations per Capita on Log Highway Gas Use per Capita"
		local noteString = "\hline \hline \end{tabular} \vspace{0.2cm} \begin{tablenotes} \begin{spacing}{0} \begin{footnotesize} \item \textbf{Notes:} The top panel of this table presents IV estimates and 95\% confidence sets robust to weak instruments of the linear relationship between log highway gas use per capita and log vehicle registrations per capita.  The bottom panel presents Anderson-Rubin chi-squared statistics along with p-values for different null hypotheses pertaining to the parameter estimated in the top panel.  The unit of observation is state-year.  The four specifications pertaining to different levels of ``Data Source Controls'' (DSCs) and control specifications are documented in Table \cref{ssec:Spec}. Confidence sets are based on standard errors clustered by state. * p\textless 0.1, ** p\textless 0.05, *** p\textless 0.01 \end{footnotesize} \end{spacing} \end{tablenotes} \end{table}"
	}
	
	local saveFile = "`tableSaveDir'"+"`tableLabel'"+".tex"
	
	local fmtDiD = `"drop(*) nocon label compress starlevels(* 0.1 ** 0.05 *** 0.01) nodepvar nomtitle nonumbers not nostar noline"'
	
	
	
	
	if("`tableLabel'"== "GasUseiv" || "`tableLabel'"=="GasUsePerCapitaiv"){
		local greekLetter = "\gamma"
		
		
		
		esttab `models' using "`saveFile'", drop(*) stats(panel1 modelNum betaPoint cs95 dsc groupFE timeFE timeTrend commonControls obs panel2 primaryARStat primaryP space1 altARStat altP space2 altARStat2 altP2 space3 altARStat3 altP3 space4 altARStat4 altP4, labels( "\hline \hline \multicolumn{5}{l}{\textit{Panel A: IV Estimates and Robust 95\% Confidence Intervals}} \\" "\multicolumn{5}{c}{Dep. Var.: Log Number of Vehicle Registrations} \\ " "$`greekLetter'$" "95\% Conf. Set" "\makecell[l]{Data Source \nextline Controls}" "State Fixed Effects" "Year Fixed Effects" "\makecell[l]{State-Specific Linear \\ Time Trends}" "Common Controls" "Number of Obs." "\hline \hline \multicolumn{5}{l}{\textit{Panel B: Anderson-Robin Test Statistics and P-Values}} \\" " \multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=11\$}} \\ AR Stat" "p-value" " " " \multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=9.7\$}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=\$1}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=\$-4.8}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=\$-6.1}} \\ AR Stat" "p-value")) nostar nomtitle nonumbers not nolines title("Robust IV Estimates of `titleOutcome' \label{tab:`tableLabel'}") replace postfoot("`noteString'") sub(\_ _)
		
	}
	else{
		local greekLetter = "\beta"
		
		esttab `models' using "`saveFile'", drop(*) stats(panel1 modelNum betaPoint cs95 dsc groupFE timeFE timeTrend commonControls obs panel2 primaryARStat primaryP space1 altARStat altP space2 altARStat2 altP2 space3 altARStat3 altP3 space4 altARStat4 altP4, labels("\hline \hline \multicolumn{5}{l}{\textit{Panel A: IV Estimates and Robust 95\% Confidence Intervals}} \\" "\multicolumn{5}{c}{Dep. Var.: Log Number of Vehicle Registrations} \\ " "$`greekLetter'$" "95\% Conf. Set" "\makecell[l]{Data Source \nextline Controls}" "State Fixed Effects" "Year Fixed Effects" "\makecell[l]{State-Specific Linear \\ Time Trends}" "Common Controls" "Number of Obs." "\hline \hline \multicolumn{5}{l}{\textit{Panel B: Anderson-Robin Test Statistics and P-Values}} \\" "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=1.6\$}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=1\$}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=0.3\$}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=0\$}} \\ AR Stat" "p-value" " " "\multicolumn{5}{c}{\underline{\$ H_{0}:`greekLetter'=-2.5\$}} \\ AR Stat" "p-value") fmt(%4.0fc)) nostar nomtitle nonumbers not nolines title("Robust IV Estimates of Log Registration on `titleOutcome' \label{tab:`tableLabel'}") replace postfoot("`noteString'") sub(\_ _)

		
	}
	
	
	
	
end
